Analytic structure of radiation boundary kernels for blackhole perturbations 
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Exact outer boundary conditions for gravitational perturbations of the Schwarzschild metric fea- 
ture integral convolution between a time-domain boundary kernel and each radiative mode of the 
perturbation. For both axial (Regge-Wheeler) and polar (Zerilli) perturbations, we study the 
Laplace transform of such kernels as an analytic function of (dimensionless) Laplace frequency. 
We present numerical evidence indicating that each such frequency-domain boundary kernel admits 
a "sum-of-poles" representation. Our work has been inspired by Alpert, Greengard, and Hagstrom's 
analysis of nonreflecting boundary conditions for the ordinary scalar wave equation. 
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I. INTRODUCTION 

Now 50 years old, the perturbation theory of Schwarz- 
schild blackholes remains a timely subject with fun- 
damental applications. The covariant d'Alembertian 
(or wave equation) associated with the Schwarzschild 
line-element describes scalar "perturbations." Classical 
Schwarzschild blackholes are spherically symmetric and 
static in time, and these symmetries allow for combined 
multipole and Fourier (or Laplace) decompositions. As 
a result, perturbations may be described via a denumer- 
able collection of ODE rather than the d'Alembertian 
PDE. Similar ODE describe electromagnetic "perturba- 
tions" U and gravitational perturbations (small genuine 
fluctuations in the background geometry) 0,01, although 
their derivation is more complicated since the multipole 
decomposition involves either vector or tensor spherical 
harmonics. Remarkably, via the technique of "despin- 
ning" based on the 3 operator £|, such ODE can be re- 
lated to scalar wave equations. 

Wheeler considered the case of electromagnetic per- 
turbations in 1955 p|, showing for a given multipole that 
each of the two electromagnetic polarization states are 
described by a copy of a single ODE. Regge and Wheeler 
then derived a similar ODE describing odd-parity (or ax- 
ial) gravitational perturbations in 1957 0, and Zerilli in- 
troduced an ODE describing even-parity (or polar) grav- 
itational perturbations in 1970 Q- In the 1970s Chan- 
drasekhar and Detweiler demonstrated that the Zerilli 
equation can be derived from the Regge-Wheeler equa- 
tion, although the derivation involves differential opera- 
tions (see |5| and references therein). In their treatment 
of intertwining operators, Anderson and Price clarified 
the relationship between solutions to the Regge-Wheeler 
and Zerilli ODE. Application of a first-order differen- 
tial operator transforms smooth solutions of one equation 
into solutions of the other. 

Schwarzschild perturbation theory has played a central 
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role in several modern areas of classical and quantum 
gravity. Although the following is by no means an ex- 
haustive list, we mention four salient applications: the 
"close-limit" approximation for blackhole collisions, a 
time-domain approach to the radiation reaction prob- 
lem, the asymptotic form of high frequency quasinormal 
modes, and quantum uncertainty in blackhole horizons. 
Price and Pullin, assuming that the colliding blackholes 
are initially cloaked in a common horizon, have used 
first-order perturbation theory and the Zerilli equation 
to compute the energy radiated away by gravitational 
waves JjJ . They initially applied their close-limit approx- 
imation to time-symmetric initial data. More general 
data was subsequently considered 0, and the accuracy 
of the approximation studied via second-order perturba- 
tion theory [jj 0] . Lousto has studied the binary radia- 
tion reaction problem in the extreme mass ratio limit via 
time-domain calculations |11| . His approach relies on nu- 
merical simulation of Schwarzschild perturbations, and 
he has developed a fourth-order convergent algorithm 
for such simulations |l2j| . Interest in the quasinormal 
mode spectrum of classical Schwarzschild blackholes has 
been renewed by possible connections with the Barbero- 
Immirzi parameter in loop quantum gravity |l3l Il4l Il5| . 
While these issues are beyond us, we note that they have 
focused attention on the asymptotic form of high fre- 
quency quasinormal modes |lo . Il7| (a purely classical is- 
sue). Most recently, York and Schmekel have considered 
a truncated superspace of blackhole fluctuations, using 
path integral quantization to derive a Schrodinger equa- 
tion and estimate the quantum uncertainty in the horizon 
[To|. York and Schmekel's analysis constitutes and im- 
proved derivation of results already found by York via 
semiclassical arguments. 

This paper addresses the mathematical issue of exact 
radiation outer boundary conditions ROBC for Schwarz- 
schild blackhole perturbations. Our work is based on 
an approach [2(j developed by Alpert, Greengard, and 
Hagstrom (AGH) for nonreflecting boundary conditions 
and time-domain wave propagation on flat spacetime. 
Beyond being of theoretical interest, such boundary con- 
ditions are relevant for numerical simulations. Indeed, 
long-time simulations require some form of domain re- 
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duction, that is specification of appropriate boundary 
conditions at the "edge" of the computational domain. 
Exact outer boundary conditions for gravitational per- 
turbations of Schwarzschild blackholes feature integral 
convolution between a time-domain boundary kernel and 
each angular mode of the perturbation. For both axial 
(Regge- Wheeler) and polar (Zerilli) perturbations, wc 
study the Laplace transform of the such kernels as an 
analytic function of (dimcnsionlcss) Laplace frequency 
a. We present numerical evidence indicating that each 
such gravitational boundary kernel admits a "sum-of- 
poles" representation. These representations are similar 
to those considered by AGH for wave propagation on fiat 
3+1 and 2+1 spacetimes. We hope to interest analysts 
in our conjectured sum-of-poles representation and the 
theorems we believe are lurking behind it. 

We have considered such boundary conditions in de- 
tail before pll |22| (hereafter referred to as Papers I and 
II). However, this paper goes beyond these references in 
the following ways. First, we consider the Zerilli equa- 
tion, not considered at all in Papers I and II. Second, 
we focus here on the analytic structure of gravitational 
kernels (both Regge- Wheeler and Zerilli cases), whereas 
Paper I almost exclusively considered kernels for scalar 
wave propagation. Although the theory and methods 
in I and II were also spelled out for the gravitational 
(spin 2) Regge-Wheeler equation, we consider these re- 
sults in more detail here. We also discuss the asymptotic 
agreement between Schwarzschild ROBC and flatspace 
nonreflecting boundary conditions, remarking on some 
fine points unmentioned in our earlier work. Part of this 
paper roughly parallels Section 3 from Paper I, which 
documented several numerical tests of the sum-of-poles 
representation for scalar (spin 0) kernels. Here we run 
through those same tests for the gravitational kernels, 
and also carry out another test based on the Argument 
Principle. A longer paper on numerical implementation 
of these ideas, complete with extensive numerical tables, 
is forthcoming |23j, and in part this paper is meant to 
lay more groundwork for that longer work. 



where F(p) = 1 — P • As is well known (see, for ex- 
ample, the review |2J|), radiative perturbations of the 
Schwarzschild metric arc ultimately described by wave 
equations for angular modes ^i m (r, p). From now on, 
we drop the "azimuthal" index m and write ^(r, p), 
since all equations depend only on the "orbital" index 
I. In terms of the Regge-Wheeler tortoise coordinate 
p* = p + log(p — 1), the perturbation obeys the wave 
equation 



dr 2 dpi 



V(p)$ e = 0. 



(2) 



For axial perturbations V(p) is the Regge-Wheeler po- 
tential 



V RW ( P ) =11- 



1(1 +1) 1 - ;f 



(3) 



where j = 2. Scalar and electromagnetic perturbations 
correspond to j — and 1, respectively. 

Eq. J5J also describes polar perturbations, although in 
this case V(p) is the Zerilli potential 



V z (p) = 



(4) 



8n 2 (n+ l)p 3 + I2n 2 p 2 + 18np + 9 



p 3 (2np + 3) 2 



where n = \ l){l+2). Throughout this paper we con- 
sider the same objects for the Regge-Wheeler and Zerilli 
cases, and we shall make use of the following convention. 
An object without a superscript letter, say V(p), will re- 
fer to the generic object, and could correspond to either 
of the two cases (actually four cases, since the Regge- 
Wheeler scenario is three cases in itself). Superscript 
letters will denote the specific cases. For example, we 
have V RW (p) and V z (p) as above. We will also sometime 
use a superscript F to denote corresponding flatspace ob- 
jects, for example 



V F {p) 



£(£ + !) 



(5) 



II. PRELIMINARIES 



A. Wave equations for gravitational perturbations 



as the flatspace potential. 

Formal Laplace transformation of the Regge-Wheeler 
wave equation (J2J yields a second-order ODE 



In terms of dimensionless time r and radius p > 1 , the 
standard Schwarzschild line-element reads 1 



ds 2 = — Fdr 2 + F~ 1 dp 2 + p 



), (1) 



1 With M the mass parameter of the solution, standard physical 
coordinates are then t = 2Mt and r = 2Mp. We also consider 
a dimensionless Laplace frequency a related to the physical fre- 
quency s via a = 2Ms. 



L RW ^ e = 0, L 



RW 



dp 



2 V RW (p) - er 2 , (6) 



where a is Laplace frequency. This equation is — apart 
from a transformation on the dependent variable — a spe- 
cial case of the confluent Heun equation pM |2^ | . There- 
fore, implementation of ROBC for the Regge-Wheeler 
equation involves confluent Heun functions, whereas the 
flatspace implementation of AGH [2(j involves Bessel 
functions (closely related to confluent hypergeometric 
functions). 
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Likewise, we may consider formal Laplace transforma- 
tion of the of the Zcrilli wave equation (J5J), 



L Z V ( = 0, 



dp* 



(7) 



Solutions of © with j — 2 are related to solutions of Q 
and vice versa by the intertwining relations 



D + L RW = L Z D + , D-L Z = L RW D-, (8) 



with 



D± = — ± 
dp* 



}n(n + l) + 



3(p-l) 



p 2 (3 + 2np) 



(9) 



We have mentioned that Papers I and II examined the 
exact ROBC for the Regge-Wheeler equation, although 
mostly focusing on the j = case. A natural question 
then is whether or not these boundary conditions can be 
easily carried over to the Zerilli case via use of intertwin- 
ing relations. It would seem the answer is "no," and we 
have been unable to make direct use of our previous work 
on ROBC via the intertwining relations. Although we be- 
lieve the issue may merit further study, in this paper we 
develop and describe Zerilli ROBC from scratch. 



B. Radiation outer boundary conditions 

Let us assume initial data of compact support, and 
that the radial location p B of the outer boundary B lies 
beyond the support of the data. Then the exact nonlocal 
ROBC is the following differential-integral identity: pH 



8t <9p* 
F(p B ) 

PB 



P=PB 



u e (T-T';p B )y e (T',p B )dT'. (10) 



The ROBC equates an outgoing characteristic derivative 
of the field with an integral convolution of the field his- 
tory. Although the conditions IjlUfl are equally valid for 
the Regge-Wheeler and Zerilli cases, we must consider 



two separate kernels: 



,RW 



(t;pb) and uf(r\p B ). As 



a function of complex Laplace frequency a, Paper I has 
considered the Laplace transform 6j RW (a; p B ) of the inte- 
gral kernel ui RW (r; p B ), in particular arguing that it ad- 
mits a "sum-of-poles" representation quite similar to the 
representation of frequency-domain kernels for flatspace 
wave propagation in 2+1 dimensions. Such representa- 
tions involve both a finite pole sum as well as a continuous 
sector J2jj. The analysis in Paper I mostly concentrated 
on the u> RW (cr; p B ) relevant for scalar wave propagation 
(that is for j = 0), and did not consider the Zerilli equa- 
tion at all. Below we consider the analytic structure of 
both wf w (cr; p B ) for j = 2 and Cbf{cr\ p B ) in some detail. 

Paper I has developed numerical methods for eval- 
uating LUg(a; p B ) along the axis of imaginary Laplace 
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FIG. 1: Reuj? w (iy, 15) and lmui^ w {iy; 15) profiles for the 
j — 2 frequency domain kernel (2>2 (iy; 15). As indicated, 
t = 2 and p B = 15. 



frequency, precisely the contour over which the inverse 
Laplace transform is taken to obtain the time-domain 
kernel. We further touch upon these methods below, but 
mention here that they rely on stable numerical integra- 
tion over various paths in both the complex p and com- 
plex z = op planes. All of our work in this paper is based 
upon these methods, and we have used them to plot in 
FIG. □ the j = 2 Regge-Wheeler profiles Reu}^ w (iy; 15) 
and Imuj2 W {iy, 15) for real y. Although different, the 
corresponding Zcrilli profiles would be indistinguishable 
to the eye, were they also plotted in FIG. ^ 



C. Kernel compression 

With the ability to generate such numerical profiles 
for exact frequency-domain kernels, we have employed 
the technique of kernel compression in order to construct 
highly accurate numerical kernels which allow for effi- 
cient evaluation of the convolution appearing in (|10() . In- 
troduced by AGH |2Cj and described further in both [33j 
and Paper II, compression is vital both for high-£ kernels 
as well as low-^ kernels which are dominated by costly 
continuous sectors (such sectors are further described be- 
low). The technique produces a rational function, 



fe=i 



7i,fc(Ps) 



(11) 



which approximates Qi(o;p B ) and is in fact a sum of d 
simple poles. The pole locations Pe,k(pB) and strengths 
H,k{pB) — output from the compression algorithm — lie 
in the lefthalf plane. The approximation is rigged to 
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k 


Re^f (15) 




Im/3«r(15) 


1 


-3.75616176922E- 


-01 





2 


-2.52285897920E- 


-01 





3 


-1.71458781191E- 


-01 





4 


-1.16562490243E- 


-01 





5 


-7.64331990276E- 


-02 





6 


-4.68091057981E- 


-02 





7 


-2.63730137927E- 


-02 





8 


-1.25652994567E- 


-02 





9 


-9.47795178947E- 


-02 


5.99312024947E- 



k 


Re/3 2 z fe (15) 




Im/3 2 z fe (15) 


1 


-3.70827103177E- 


-01 





2 


-2.48916278532E- 


-01 





3 


-1.69097983283E- 


-01 





4 


-1.14894340803E- 


-01 





5 


-7.53169595130E- 


-02 





6 


-4.61252633339E- 


-02 





7 


-2.59880681276E- 


-02 





8 


-1.23819599759E- 


-02 





9 


-9.34065839850E- 


-02 


5.89802789971E- 



A: 


Re 72 T(15) 




Im 72 T(15) 


1 


-9.42815440764E 


-06 





2 


-3.66046310052E- 


-04 





3 


-3.74027383589E- 


-03 





4 


-8.72734265927E- 


-03 





5 


-1.47189136342E 


-03 





6 


-5.01356988668E 


-05 





7 


-9.73423621068E 


-07 





8 


-7.28807025058E 


-09 





9 


-8.94836172991E 


-02 


6.20643548937E-02 



k 


Re 72 z fe (15) 




Im7a^(15) 


1 


-8.94981130535E- 


-06 





2 


-3.47903587670E- 


-04 





3 


-3.56001528614E- 


-03 





4 


-8.35011768248E- 


-03 





5 


-1.41306545024E- 


-03 





6 


-4.81777344134E- 


-05 





7 


-9.35693850095E- 


-07 





8 


-7.00652022551E- 


-09 





9 


-8.70154844197E- 


-02 


6.01803999946E- 



TABLE I: Compressed j = 2 Regge- Wheeler kernel for I = 2, 
pB = 15, £ = 10 -10 . There are d = 10 poles and strengths, 
and complex conjugation of the ninth entries gives the tenth 
entries. Zeros correspond to outputs from the compression 
algorithm which are less than 10 -30 in absolute value. 



TABLE II: Compressed Zerilli kernel for^ = 2, p B = 15, 
e = 10 -10 . There are d = 10 poles and strengths, and complex 
conjugation of the ninth entries gives the tenth entries. Zeros 
correspond to outputs from the compression algorithm which 
are less than 10 -30 in absolute value. 



satisfy 



D. Quasinormal ringing and decay tail 



\ui(iy;p B )\ 



< e, 



(12) 



where e is a chosen numerical tolerance. Theoretically, 
this bound on the relative supremum error in the fre- 
quency domain ensures a long-time bound on the rela- 
tive convolution error associated with (|10|) . as discussed 
in |20l | and Paper II. This convolution error arises when 
using the approximate time-domain kernel 



€t(r; Pb) = ^2 li,k{PB) exp [Pi,k{p B )T] 



(13) 



k=l 



in place of the true kernel uji{t\Pb)- Since ^(t;pb) is 
a sum of exponentials, due to recursive identities an ap- 
proximation of the convolution IjlUII based on Ci( t 'tPb) 
is not memory intensive. We list representative com- 
pressed kernels for Regge- Wheeler and Zerilli kernels in 
Tables HI and UTI The two kernels are strikingly similar, al- 
though they differ by a relative supremum error of about 
3.03 x 10 -2 , well below their stated e = 10~ 10 tolerances. 
Extensive numerical tables of compressed kernels are be- 
ing prepared and will appear in 23] . 



Providing little detail, we now carry out a simple 
experiment meant only to indicate that the described 
ROBC work well for long-time simulations. More care- 
ful experiments were considered in Paper II (for Regge- 
Wheeler cases only) and will be considered in |2^. In 
terms of retarded time /i = r — p* and the pulse function 



s(m) 



[ M ( M + 4)] 4 /256, 
0, 



for - 4 < a* < 0, 



otherwise, 
we construct the t = 2 wave packet 

*2(o,p) = 9(-p*) 



(14) 



(15) 



as initial data for the Zerilli equation. The initial packet 
is then of unit height and compactly supported on < 
p* < 4. To complete the data, we assume that 



<9t 



dp* 



= 0, 



(16) 



T = 



so that the pulse starts as essentially outgoing. We place 
the inner boundary at /?* = —175, and adopt (|10fl as the 
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nates, although for the flatspace case r, t, and s would 
be more standard notations. Their boundary condition 
for a flatspace order-£ multipole is 



cm 

dp 



99 t 
dr 

i r Ft 

PB Jo 



r';p B )^(r',p B )dr', (17) 



comparable with (|10|l . The exact TDRK uif(T;pB) ap- 
pearing in (|17|) is now also a time-domain nonreflecting 
kernel (TDNK) , and it is the inverse Laplace transform of 
an exact frequency-domain nonreflecti ng k ernel (FDNK) 
admitting the following representation: [20j 



FIG. 2: Quasinormal ringing and decay tail. 



boundary condition at p = ps = 15, with an approxima- 
tion of the exact ROBC based on the compressed kernel 
listed in the Table |H| 

We evolve the data until r = 300, using the MacCor- 
mack predictor-corrector algorithm. 2 Throughout the 
evolution, we record the value *& 2 (t, 3.25) of the field, 
that is to say, we record the history of the field in time 
at the fixed location p — 3.25 (actually at the grid point 
nearest this location). Notice that the total run time 
is not long enough for the history if? 2 (r, 3.25) to be in- 
fluenced by reflection off of the inner boundary. The 
absolute value of the history is depicted as a linear-log 
plot in FIG. |21 It exhibits quasinormal ringing |24j until 
about r = 150, and afterward a Price tail with the field 
decaying as i~( 2f + 3 ) = t~ 7 |24|. This late-time behavior 
stems from backscatter of the outgoing packet off of the 
long-range potential V z {p). The dashed curve 3000i~ 7 
has been eyeballed to fit the late-time decay tail. 



III. SUM-OF-POLES REPRESENTATION 

A. Flatspace FDNK and outgoing solution 

AGH have considered nonreflecting boundary condi- 
tions (NRBC) for both 3+1 and 2+1 flatspace wave prop- 
agation, thoroughly treating both theoretical description 
and numerical approximation of NRBC for both scenar- 
ios |2(|. For later comparison with the two gravitational 
scenarios, let us briefly recall the principal theoretical 
aspects of their work for the 3+1 scenario. To facili- 
tate the comparison, we will use the same letters p, r, 
and a used for our dimensionless Schwarzschild coordi- 



uf{o-,p B ) 



k=l 



be,k/pi 



(18) 



where the bg^ are the zeros of the classical MacDonald 
function -fQ +1 / 2 (z), a modified cylindrical Bessel func- 
tion. Here the Bessel order is a half-integer, since we 
are considering the radial wave equation stemming from 
ordinary wave propagation on 3+1 flat spacetime, and 
these functions have the form 1331 



K V2 (z) 
K 6/3 (z) 
K 7/2 (z) 

K e+ i/ 2 (z) 



2z 

7T 



2z 

7T 



2z 

7T 



1TZ 
IP 



1 

1 + - 

Z 



1 3 3 

H h 

z z z 



6 15 15 

H h + ^7 

z z z z a 



■i)V[H) e -. 

z az J z 



(19) 



The function K i+1 / 2 (z) has t simple zeros : k — 

1, • • • , £}. When scaled by order, these zeros (£ + 
l/2)~ 1 bi i k are known to accumulate on a fixed transcen- 
dental curve C in the lefthalf plane [2(], 13 (^ ne numeri- 
cal methods developed in Paper I take advantage of this 
fact). Although this accumulation is asymptotic with 
large order I + 1/2 — > oo, FIG. shows that the agree- 
ment holds even for the lowest I, at least to the eye. The 
curve C shown in FIG. has parametric form pMl^ 



z(A) = - VA 2 - A tanh A + iyA coth A - A 2 



(20) 



2 A scheme in consistent conservation form '27] when applied to a 
hyperbolic conservation law. 



for A in the domain [0, Ao] with Ao — 1.1997 such that 
tanhAo = I/Aq. In terms of the "normalized-at-infinity" 
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FIG. 3: Scaled zeros of MacDonald functions. Here we 
plot scaled zeros (£ + l/2) -1 6^fc for £ = 1,2, 3, 4. The cross is 
the scaled zero of Ki/2(z), the diamonds are the scaled zeros 
of K 3 / 2 ( z ), the circles are the scaled zeros of K s / 2 (z), and the 
stars are the scaled zeros of K7/2 (z). 



outgoing solution 3 



W e (z) 



exp(z)K e+1 / 2 (z), 



we have 



W' p {ap B ) 
Wi{(xp B ) 



(21) 



(22) 



as another expression for the flatspace frequency-domain 
kernel 20]. In passing, we remark that although the ex- 
act FDNK (|18fl is already a rational function, the tech- 
nique of kernel compression still proves useful for high-£ 
FDNK, since for a given tolerance e compression of <|18[) 
yields a numerical kernel with far fewer poles |20| . 



B. Gravitational FDRK and outgoing solution 

Let us first consider the relationship between a gravi- 
tational FDRK uig(a; ps) described earlier and solutions 



3 In the introduction of Paper II, the correspondence between 
Wi(z) and -fQ+i/^C- 2 ) is off by a factor of tt/2. 



to the formal Laplace transform 



d 2 fg 
dpi 



V{p)^ t 



(23) 



of the generic equation @. Here we mostly just collect 
relevant formulas. These formulas are derived in the first 
section of Paper I, and the derivation presented there 
goes through for Regge- Wheeler case (j — 0,1,2 con- 
sidered there) as well as the Zerilli case (not considered 
there). Whence the formulas we consider now are valid 
for all cases. 

We "peel off" the exponential behavior of the field by 
setting it £ — exp(— crp*)$^, whereupon finding 



"dy" 



-2a 



1 

P 



1-2(7 



dp 



P 2 V(p) 



*# = 



p-1 J dp (p-l)a 

(24) 

as the ODE satisfied by When numerically integrat- 
ing (|24|1 in various contexts, we have found it useful to 
work instead with z = ap, re-expressing the equation as 
follows: 



d 2 f> e 
dz 2 



dz 



a 2 {z-af 



(25) 

Let We(z;<j) denote the outgoing solution to l|25|l . This 
solution obeys |28| 



(26) 



k=0 



as z — > 00, and we describe it as "normalized at infinity." 
Appendix B considers the recursion relations defining the 
dfr W (a) and <if (a). As described in Paper I, the flatspace 
solution Wt(z) in JHJ) is formally W e (z) = W e RW (z; 0) in 
terms of the Regge- Wheeler case W e RW (z; a). 

In our notation Wi(ap;a) is the outgoing solution to 
(I24|l . and Paper I expresses the FDRK in terms of it as 



iot\a\p B ) 



opb 



W[{ap B \o) 
We{apB]cr) ' 



(27) 



where the prime denotes differentiation in the first slot 
of Wi(z;a). Since we have peeled off the exponential 
factor exp(— up*) above and now work with an outgoing 
solution Wg (ap, a) which is normalized at infinity, the 
FDRK (|27|l is more apt to have a well-defined inverse 
Laplace transform. 

Paper I has described a collection of numerical meth- 
ods which allow us to (i) evaluate Wg (&Pb \ o) for a in 
the lefthalf plane and (ii) evaluate togiiy] Pb) for real y. 
We remark that Leaver has analytically represented a 
solution, say the outgoing one Wj RW (apB', c) exp(— ap*), 
to the frequency-domain Regge- Wheeler equation as an 
infinite series in Coulomb wavefunctions, where the ex- 
pansion coefficients obey a three-term recursion relation. 
In fact, such series expansions hold more generally for 
the generalized spheroidal wave equation (essentially the 
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confluent Heun equation) [2!j ■ The methods described in 
Paper I are not based on the appropriate Leaver series, 
rather they rely on direct integration of (|25p. As such 
they can and have now been carried over to the Zerilli 
case. We note that the Zerilli equation is not directly 
related to the confluent Heun equation; whence it is not 
immediately evident how to evaluate Wf (opB', &) via a 
Leaver series. 

Roughly, our method for evaluating Wt{ops]Cf) is as 
follows. For a fixed frequency o, initial data at a large 
radius are obtained for i|25|) using the asymptotic expan- 
sion (|26l) . Then l|25l) is integrated over a suitable path in 
the complex z — op plane to a terminal point zb = ops- 
To convey the basic idea behind our numerical evaluation 
of the FDRK oje(iy;pB) itself, we introduce 



wt(z; a) = z 



(28) 



Then the FDRK is (2>g(o; pb) — we(op B ', c), and we have 
from H25fl that wi(z;o) obeys the first-order nonlinear 
equation 

dwi wf ( 2 Ma\ z 3 V(z/o) 

- 1 h -2 1 we 2? = 0. 

az z \ z z — oj o z (z — o) z 

(29) 

We calculate values uii(\y;pB) = u!i(iyps',iy) via direct 
numerical integration of 1|29|) . To achieve stability and 
high accuracy, the methods associated with both (i) and 
(ii) evaluations require integration over nontrivial paths 
in the z-plane. Moreover, for technical reasons the inte- 
gration employed for kernel evaluation (ii) is sometimes 
carried out in the complex p-plane rather than complex 
z-plane. Finally, we remark that to produce the origin 
value LUi(0;ps) of the kernel, we do not integrate (|2*9")l . 
Rather, we make use of an exact series expression, one 
given for the Regge- Wheeler cases in Paper I and for the 
Zerilli case in Appendix A. 



depends on the value of ps - 4 In principle the integer Ng 
also depends pb, but turns out to be constant over sizable 
regions of the relevant parameter space (more comments 
on this point below). The A;th pole strength and cut 
profile are given respectively by 



fe(Xi Pb) = lmw e (xe 17T ; Pb) , 
(31) 

with x > and the prime here standing for d/dps dif- 
ferentiation. 

We will argue that the representation (|3(J[1 is also valid 
for both j = 2 Regge- Wheeler and Zerilli gravitational 
cases, and this conjecture is the main result of our paper. 
Note that lJ5U|l is not really a "sum of poles." Indeed, re- 
call that in the sense of complex analysis a pole is an iso- 
lated singularity. Strictly speaking then, the cut integral 
in (|3U[) does not correspond to a "continuous distribution 
of poles" (an oxymoron). Nevertheless, we shall continue 
to describe the representation (|30l) as a "sum of poles." 

We stress that (|30|l is a representation for a bound- 
ary integral kernel. Its Laplace transform, the TDRK 
loi>(t] pb), lives on the history of the spatial boundary B. 
Indeed, (|10|) makes no reference to the details of the ini- 
tial data and is certainly not a spatial convolution over 
initial data. Moreover, the pole locations ot,k{pB) are not 
quasinormal modes. An infinite number of quasinormal 
modes belong to each £ value, and these characteristic 
frequencies do not depend on any particular choice of 
outer boundary radius ps- Quasinormal modes are as- 
sociated with a boundary value problem specifying that 
^ i is downgoing at the horizon and outgoing at infinity. 
The locations o~e,k(pB) are finite in number, and they do 
depend on ps- They can be associated with a boundary 
value problem specifying that ^ g is outgoing at infinity 
and vanishes at pB- b Likewise, the cut integral appearing 
in H3()|) is not the branch-cut contribution to the usual 
Green's function studied in the quasinormal mode prob- 
lem El. 



C. Gravitational sum— of-poles representation 

For the case of j = scalar perturbations Paper I 
has documented compelling numerical evidence indicat- 
ing that the FDRK (|27|l admits an explicit "sum-of- 
poles" representation, 

ui(o\ pb) = > — r - - / dx , 

^ er- °-f„k(PB) 7T J er + x 

in terms of complex frequency a. When we discussed 
compressed kernels before, we introduced approximate 
pole locations Pe t k(pB) and strengths "fe^ipB)- We now 
consider Ni (an integer) physical pole locations og^iPB) 
and physical pole strengths ae^k(pB), all of which lie in 
the lefthalf plane. Also appearing in l|3U|) is a cut profile 
fi{XiPB), and like the pole locations and strengths it 



IV. NUMERICAL STUDY 

In this section we both provide further qualitative de- 
scription of the key representation i|30|) and justify it nu- 
merically. To argue that the representation (|30|l is valid 
for both j = 2 Regge- Wheeler and Zerilli gravitational 
cases, we offer nearly the same evidence as that offered 
in Paper I for the case of scalar perturbations. However, 
here we also consider one extra numerical experiment 



4 It may or may not be the case that an approximate location 
Pl.k(PB), for example, approximates a physical one ut,k{pB)- It 
is the compressed kernel &(iy] ps) in whole which approximates 
the physical FDRK iOe(iy; pb) uniformly in y 6 K. 

5 As such, the locations k(PB) arc analogous to the "fiatspace 
quasinormal modes" considered in 131,1 , a misleading terminology 
for our paper. 
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based on the Argument Principle. Our main numerical 
justification is to compare values u>e{iy, pb) of the kernel 
obtain via two independent approaches. These are the 
following: 

• Direct integration of l|29|) as alluded to above [a 
process which makes no use whatsoever of the con- 
jecture representation J3DJ]. 

• Approximation of the sum-of-poles representation 
itself (by this we do not mean kernel compression) . 

In the second approach, we "build" the kernel as 



ui{w,Pb) 



<*i,k{PB) + Q~ai,k{pB) 

f^Ay — o-£,k(pB) - 5o- Lk (p B ) 

_ l f x — Mx;pb) + 6Mx;pb) 
77 h^r, w + x 



(32) 



&X- 



where the S terms represent numerical errors and the in- 
tegral over the chosen window [xmiiu Xmax] must be han- 
dled via numerical quadrature. We have used Simpson's 
rule, respectively with 2048, 2048, 1024, 2048, 1024, 2048, 
1024, 2048, 512 subintervals for I = 2, 3, 4, 5, 6, 7, 8, 
9, 10 (odd I require more subintervals). We stress that 
this latter approach to evaluation, quite unlike the first 
approach, requires that we numerically compute all pole 
locations and strengths and well as the cut profile. To 
locate poles, zeros in a of Wi{aps] c), we have used the 
secant algorithm. The d/dps derivatives of the <Jt, k {pB) 
needed to compute the strengths a.i,k{.PB) are obtained 
by first building a high-order interpolating polynomial 
Ti^{1/pb) for each <T£^(pb) based on Chebyshev nodes 
in 1/ Pb- Derivatives are then found via differentiation of 
the Chebyshev polynomial. This procedure is described 
in more detail in Paper I where it was used for the j = 
Regge- Wheeler case. 

We remark that this direct approximation (|32|l to the 
sum-of-poles representation (|30|l is certainly not a com- 
pressed kernel. Indeed, as mentioned, this brute-force 
approximation of (|30|l typically requires thousands of 
poles (stemming from the numerical quadrature of the 
cut integral) to achieve the same error tolerance achieved 
by a compressed kernel comprised of ten or so poles. 

Besides making the comparison outlined in the last 
paragraph, we also wish to provide further qualitative 
description of the sum-of-poles representation l|30|l . We 
describe both the poles and the cut profile in more de- 
tail, and also compare the representation to the strikingly 
similar representation (|18f) of the FDNK ui£ (er; ps) for a 
flatspace order-i? multipole. We have found that both 
gravitational kernels, Cbf w (a; pb) and u>f(a; pb), indeed 
agree with ui[(o~;pb) in the ps — > oo limit, although 
here we will focus on the Zerilli case. One certainly ex- 
pects such agreement, since the Schwarzschild solution is 
asymptotically flat. However, the nature of this asymp- 
totic agreement is rather interesting, and we point out 
some subtleties not mentioned in Paper I. 




FIG. 4: Zeros {cr^ fc (15) 



Reo 



1 < k < 4} in frequency of 



Wg (<t15;o-). The contour lines are of log 10 | W£ (crl5; a)\, with 
the logarithm distributing contour lines more evenly. 



0.15 




Reo 



FIG. 5: Zeros in frequency of IF;f (crl5; a). The plot here is a 
blow-up of the one shown in FIG. 2] 



A. Pole locations and strengths 

For the most part, this subsection gives what amounts 
to a qualitative description. Using our method for eval- 
uating the outgoing solution Wi(upB]o-), we may plot 
the modulus | Wraps', &)\ in order to suggest rough val- 
ues for zeros (that is, roots) of the function. Provided 
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such a zero is simple, it will correspond to a pole appear- 
ing in the representation (|30|l . Over the parameter space 
Pb > 15 and £ = 2, 3, • • • 10, and for both Zerilli and j = 2 
Regge- Wheeler cases, we have found that the number Ng 
of zeros <Je,k(pB) is as follows: N 2 — 2, N 3 = 4 = N 4 , 
N 5 = 6 = N 6 , N 7 = 8 = iV 8) N 9 = 10 = iVm. For a fixed 
choice of I and ps, these zeros form a crescent pattern 
in the lefthalf cr-plane. 

For example, FIG.0]depicts the modulus |W^(<rl5; tr)| 
in the indicated region of the lefthalf a plane. Four 
zeros appear to be evident in the figure, and in order 
to further explore whether or not they are indeed ze- 
ros, we appeal to the Argument Principle. We focus on 
the two upper locations shown closer up in FIG. Let 
h(a) — Wf(al5; a) represent our numerically computed 
function. 6 We numerically compute 

——7 / T7-T dcr - 1-999732 - i4.185506 x 1(T 5 (33) 

27T1 /square M ") 

over the square (running in the counterclockwise sense) 
shown in FIG. [5] On each side of the square we have in- 
troduced 1024 subintervals, and used the trapezoid rule. 
Since we are unable to numerically evaluate h'(cr) di- 
rectly, we approximate this derivative using difference 
quotients (this requires two extra function evaluations 
beyond the corners). Using a sequence of discretizations, 
we have confirmed that the integral convergences to 2 at 
a second-order rate, suggesting that the square indeed 
encloses two simple zeros. As a second example, con- 
sider the modulus | W 8 z (ct25; a)\ plotted in FIG. over 
the indicated region (four other zeros, conjugate to those 
shown, are located in the third quadrant) . For the analo- 
gous line-integral over the square shown (and using 2048 
cells on each side for the trapezoidal integration) , we find 
a value of 3.999975 - i5.169761 x 10~ 5 , suggesting four 
simple zeros. 

As with the scalar case, we find for the Zerilli and j = 2 
Regge- Wheeler cases that the zeros <Ji,k{pB,) m frequency 
of We(apB',cr) behave asymptotically as 

V£,k(PB) ~ h,k/PB (34) 

in the ps — > 00 limit. However, let us offer several impor- 
tant observations in order to sharpen the precise nature of 
this asymptotic agreement. First, in the opposite limit as 
ps — > 1 + , we do not believe Ni remains constant. Indeed, 
we expect the phenomenon of "zero pair creation" in this 



Due to small errors in the asymptotic expansion 1261 used to 
generate initial data for an evaluation based on integrating 1251 . 
h(o) will actually represent the product of (crl5; cr) and an 
analytic function of a which slowly varies over the region of in- 
terest. However, Paper I showed that the choice of integration 
path in the z = op plane results in exponential suppression of 
the second solution to 1251 . and absolute differences in the zero 
locations of h(a) and those of W^(ol5;o~) are of size 10 -13 in 
modulus. 




FIG. 6: Zeros in frequency of Wg(a25; a). The contour lines 
are of log 10 |W g z (cr25; cr)|, with the logarithm distributing con- 
tour lines more evenly. 



limit, as described in Paper I for the scalar j — Regge- 
Wheeler FDRK. Although we believe that the sum-of- 
poles representation remains valid for ps close (but not 
equal) to unity, description of the FDRK and implemen- 
tation of ROBC both become more difficult in this limit. 
For these reasons, we have required ps > 15. Second, for 
odd £ there is an "extra" zero. That is to say, for odd 
£ the number Ni = £ + 1 of zeros <Je t k(pB) is greater by 
one than the number i of MacDonald zeros (provided ps 
is large enough, otherwise Ni could be a larger integer 
still). We will argue that this curious feature is not at 
odds with the asymptotic result (|34fl above. Let us focus 
on the Zerilli case, with the understanding that similar 
statements apply to the j = 2 Regge- Wheeler case. 

For even I — 2,4,6,8,10 and ps > 15, the num- 
ber of zeros crf k (pB) is Ni — £, the same as the 
number of bi^/pB- As an illustration, FIG. \7\ depicts 
the six zeros cr^ fc (ps) of W§ '(apB^cr) as diamonds for 
ps = 15, 16, • • • ,30. We also plot the corresponding 
MacDonald-Bessel zeros be,k/pB as crosses. The col- 
lection CT|? fc (15) of zeros is the outermost crescent of di- 
amonds, while the collection <7j? fc (30) is the innermost 
(and similarly for the crescents of Bessel crosses). The 
plot is clearly not at odds with the asymptotic formula 
above. 

As mentioned, for odd £ — 1,3, 5, 7, 9 and ps > 15, the 
number of zeros crf k {pB) is Nt = £ + 1, that is one more 
than the corresponding number of MacDonald-Bessel ze- 
ros. For odd £ there is a single MacDonald-Bessel zero 
which lies on negative real axis. As with the j = 
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FIG. 7: Zeros {erf \(j)b) : 1 < < 6} in frequency of 
W§(ppB\ cr). Diamonds represent the zeros cr6,fc(ps), while 
the crosses represent the be,k/pB- All zeros for ps = 
15, 16, • ■ • ,30 are shown. The outermost crescent of locations 
corresponds to ps = 15 and the innermost to pb = 30. 



Regge- Wheeler case, we find for odd I that two zeros of 
Wf (apB',cr) correspond to this single real MacDonald- 
Bessel zero. Moreover, as ps gets large each of these two 
zeros is asymptotic to the single MacDonald-Bessel zero. 
This phenomenon is evident in FIG.|H1 and correspond- 
ing plots for other odd t are similar. The existence of 
an "extra" zero is at first sight troubling in light of the 
expected asymptotic agreement between wf(iy;pB) and 
(iy; ps)- However, we argue below that the pole and 
cut contributions to the gravitational FDRK in tandem 
do yield the correct asymptotic agreement. 

So far we have only consider locations of k {pB) asso- 
ciated with Zerilli kernels. To the eye, both FIGS. 
and 03 would be the same had we instead plotted the 
corresponding locations ct^(pb) associated with j = 2 
Regge- Wheeler kernels. For the examples we have con- 



sidered, \crf k (pB) — a £k{PB)\ is typically on the order 
of 10~ 5 to 10~ 3 . However, from Paper I results we know 
that this difference is seven to ten orders of magnitude 
larger than our knowledge of these locations. As an ex- 
ample, we list the two conjugate locations 02,fc(15) for 
j = Regge- Wheeler, j — 2 Regge- Wheeler, and Zerilli 
cases. Respectively, these are 

-0.096885391711329 ± 10.061245961499841, 
-0.094779501145744 ±10.059927941363806, 
-0.093406539086840 ± 10.058977077340679. 

These locations are all quite similar, and they may be 
compared to the roots &2,fc/15: —(3/2 ± iv / 3/2)/15 = 
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FIG. 8: Zeros 

of Wf(crpB', cr)- Diamonds represent the zeros a$ k (pB 
while the crosses represent b^^/pB- All zeros for ps = 
15, 16, • • • ,30 are shown. The outermost crescent of locations 
corresponds to ps = 15 and the innermost to ps = 30. For 
each pB two diamonds correspond to the single cross lying on 
the negative real axis. 



-0.1 ± i0.05773502691896. In this short list the last lo- 
cation, say, clearly corresponds in some sense to the ninth 
pole location in the compressed kernel in Table [n] but 
note that the absolute difference between the value here 
and the table value is more than the stated e = 10~ 10 
tolerance for the compressed kernel. See the footnote just 
after Eq. PP|). 



B. Cut profile 

Using the type (i) numerical method from Paper I de- 
scribed above [which also returns derivative information 
for W^(apB',cr)], we have generated Zerilli cut profiles 
ff(X'i Pb) for various I and pb- FIG. depicts profiles 
for even I. To the eye the cut profiles shown match j = 2 
Regge- Wheeler profiles f? (x; Pb) corresponding to the 
same parameter choices. However, the Zerilli and Regge- 
Wheeler even profiles are indeed different, as is shown by 
an example plot in FIG. ^] Near the peak of this plot- 
ted difference, the accuracy to which we know the Zerilli 
and Regge- Wheeler profiles is some ten orders of magni- 
tude greater than the difference. Notice that these even 
cut profiles weaken as ps gets larger, and we believe that 
the cut contribution to an even-^ gravitational FDRK 
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FIG. 9: Scaled even cut profiles for Zerilli kernels. In each 
plot the leftmost profile is fi (x> Pb)/2-5 and the rightmost 

/fo(x;ps)/io.5. 
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FIG. 10: Scaled odd cut profiles for Zerilli kernels. In each 
plot the leftmost profile is //(xiPs)/3-5 and the rightmost 

/•f(x;ps)/9.5. 



"dies out" as ps gets large. Based on our earlier claim 
that for an even-£ gravitational FDRK we have the same 
number Nt = £ of zeros <Ti,k(pB) as Bessel zeros bg^l Pb-, 
and that these zeros lock on to the latter as Pb gets 
large, we conjecture that the ev&n-l kernels u)f (iy; pb) 
and &f w (iy; ps) both approach Cuf (iy; ps) uniformly in 
y as pb — > 00. 

We plot odd-^ Zerilli cut profiles in FIG. EJI Like be- 
fore with the even-^ profiles, to the eye these could be 
either Zerilli or j — 2 Regge- Wheeler profiles. However, 
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FIG. 11: Difference between Regge- Wheeler and 
Zerilli cut profiles. Here we plot example differences with 
the notation A/,( X ; 15) = jf w ( X ; 15) - /# ( X ; 15). 



as is evident in the bottom plot shown in FIG. ITU these 
cases are different. Notice that for odd I the cut pro- 
files strengthen as ps increases. We believe that such a 
strengthening profile in tandem with the two zeros clos- 
est to the real axis as a whole combine to asymptotically 
agree with the single MacDonald-Bessel zero be,o/pB lo- 
cated on the negative real axis. [For notational simplic- 
ity we have now labeled this zero by k = 0. Before it 
would have corresponded to k = \(l + 1) for odd £, if 
k = 1, • • • ,£ for the zeros be,k of Kg +1 / 2 (z).] In other 
words, for odd-£ the cut profile in effect cancels one of the 
zeros <Je,k(pB), so that we again we have agreement be- 
tween the gravitational FDRK Cof(iy;pB) and flatspace 
FDNK u>f(iy; ps) uniformly in y as ps — > 00. While we 
cannot precisely describe how this cancellation of the ex- 
tra zero occurs, we believe that it stems from the follow- 
ing two conjectures, claimed to hold as ps — * 00. First, 
the cut profile becomes sharply concentrated for x near 
-bifi/pB- Second, 



1 



h(x;PB)dx — h,o/pi 



(35) 



Preliminary numerical investigations indicate that both 
claims are in fact valid, but the issue deserves further 
study (preferably theoretical). All statements made in 
this paragraph also pertain to ilif w (iy; ps)- 



C. Numerical validation of the sum— of— pole 
representation 

Up to this point we have mainly given a qualitative de- 
scription of what we believe are the main features of the 
sum-of-poles representation (|30|) . In order to quantita- 
tively test our representation (|30|) . we compare the two 
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FIG. 12: Relative and absolute numerical errors. 
Here we plot numerical errors corresponding to \Auj§ (iy; 15) | 
and [Awf (iy; 15)|/|&# (iy; 15) | for the Zerilli case with t = 2 
and pB = 15. 



numerical approaches for evaluating either Cjfiiy; pb) or 
wf w (iy; Pb) outlined at the beginning of this section. 
FIG. 1121 depicts such a comparison for the Zerilli case 
with I — 2 and ps = 15. To make the plots in the figure, 
we have used a y-grid with 5 adaptive levels, each with 
32 grid points. The adaptive grid provides more resolu- 
tion near the origin where we expect the largest errors. 
We have then generated two separate numerical arrays 
of values a>;f (iy; 15) on the y-grid, the first obtained by 
integrating the ODE (|29|l and the second via the direct 
construction l|32|l . Both arrays of numerical values are 
obtained in double precision arithmetic. The plots de- 
pict the absolute error | Acj^(iy; 15)| and relative error 
|Aa)f (iy; 15)|/|a)f (iy; 15) | over the y-grid, and they in- 
dicate striking agreement between the two methods. We 
note that the maximum value of the absolute error in 
the top plot is 1.49 x 10~ 12 , and the maximum value of 
the relative error in the bottom plot is 7.35 x 10 -13 . We 
can now perform the same experiment over a range of 
Pb values, say ps = 15, 16, • • • , 30. For each choice of 
Pb we compute maximum values of absolute and relative 
error over the y-grid. It turns out that the same values 
above corresponding to pB = 15 are the largest errors 
encountered. 

We now perform the same experiment for each t < 10 
and for both j = 2 Regge-Wheeler and Zerilli cases. 
That is to say, for each £ we compute the maximum ab- 
solute error |Ac2^(iy; pb)\ and maximum relative error 
\At2)g(iy; pB)\/\u)e(iy', Pb)\ uniformly over the described 
y-grid and all ps = 15, 16, • • • , 30. For the Zerilli case 
we list these errors in the Table IIIII Errors for the j = 2 
Regge-Wheeler case are comparable. 
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7.35 x 10" 
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1.49 x 10" 
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3.63 x 10" 
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9.50 x 10" 
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8.63 x 10" 
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2.47 x 10" 
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1.27 x 10" 
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4.68 x 10" 
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1.51 x 10" 
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5.10 x 10" 
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3.51 x 10" 


12 


10 


9.71 x 10" 


14 


8.92 x 10" 
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TABLE III: Zerilli errors 



V. CONCLUSION 

For both flatspace 3+1 and flatspace 2+1 wave prop- 
agation and nonreflecting boundary conditions, AGH 
proved several theorems related to both the exact sum- 
of-poles representation for a FDNK (also a FDRK, but 
we have used FDNK for t he sp ecial flatspace case) and its 
numerical approximation [20j. In particular, they proved 
that a FDNK admits a rational approximation as a com- 
pressed kernel, and for a given large Bessel order v (for 
the 3+1 case v = I + 1/2) and choice of e tolerance they 
estimated the required number d of poles appearing in 
the compressed kernel. In addition to some techniques 
stemming from the Fast Multipole Method, their proofs 
rely heavily on the well-understood theory of Bessel func- 
tions. In particular, AGH made extensive use of integral 
representations, order recursion relations, and detailed 
understanding of the scaling behavior for both the poles 
and (for the 2+1 case) the cut profile. 7 

While the numerical evidenced amassed here is ex- 
tremely convincing, it does not constitute a mathematical 
proof that a gravitational FDRK admits the representa- 
tion H30(l . and we hope that our numerical investigation 
spurs the interest of analysts capable of theoretically in- 
vestigating our conjectured representation. Were lj5U|) 
established theoretically, we believe approximation the- 
orems — similar to those proved by AGH but pertaining 
to our gravitational ROBC — would follow. One might 
prove that a gravitational FDRK also admits a ratio- 
nal approximation as a compressed kernel, and determine 
the asymptotic growth of the number d of approximating 
poles as £ — > oo, e — > + . Since numerically constructed 
compressed kernels have performed spectacularly in im- 
plementations of ROBC (see Paper II), we have good rea- 



7 The 2+1 case, associated with integer-order Bessel functions, 
is analytically richer than the 3+1 case, as for the 2+1 case 
there is a continuous sector associated with the sum-of-poles 
representation for the FDNK. This sector is at least qualitatively 
similar to the one in 1301 . 
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son to believe that approximation theorems must hold- 
Such theorems are bound to involve the details of the 
underlying special functions. Unfortunately, relative to 
Bessel functions, significantly less is known about the 
special functions considered here, confluent Hcun func- 
tions for the Regge- Wheeler cases and seemingly more 
exotic functions for the Zerilli case. Indeed, we are un- 
aware of useful integral representations, and while appro- 
priate Leaver series are certainly of formal interest, they 
would not seem a good platform for carrying out the 
requisite asymptotic analysis. While our own knowledge 
of modern analysis would seem not up to such theoreti- 
cal investigation, we believe our results offer fertile new 
ground for more capable analysts. 
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APPENDIX A: ORIGIN VALUE OF THE FDRK 

For the j = 0, 1,2 Regge- Wheeler cases Paper I has 
expressed the origin value w^ w (0; ps) of the FDRK in 
terms of an infinite series in p~^ , where the expansion 
coefficients obey a two-term recursion relation (see Sec- 
tion 3.2.3 of that reference). All expansion coefficients 
are positive, and the value ujf w (0; pb) can be accurately 
approximated in terms of two partial sums. This series 
was used numerically to evaluate the kernel at the origin 
a = 0, a frequency value for which direct integration of 
(|29|l proved problematic. 

To express u>/(0;/jb) for the Zerilli case of polar per- 
turbations, we have used a similar series expression, al- 
though now the series coefficients obey a four-term rather 
than two-term recursion relation. We have obtained this 
series via the following recipe. First, we set a = in i|24|) 
and use ©, thereby reaching 



where 



d 2 $. 



1 



1 \ d$/ 



dp 2 V /) dp 

"8n 2 (n + l)p 3 + 12n 2 + 18np + 9 



p 2 {p~ l)(2n/7+3) 2 
Next, plugging the expansion 



O-kP 



-(e+k) 



$>e = 0.(A1) 



(A2) 



k=0 



into i|Al|l . assuming ao = 1, and balancing terms, we find 
the recursion relation 



/3 fc 
7fe 



= 4n 2 [(£ + k + 3)(£ + k + 4)-2(n + l)}, 

= in(£ + k + 2)[(3 - n){£ + k + 3) — n] — I2n 2 , 

= (£ + k + 1)[(9 - 12n)(£ + k + 2) - 12n] - 18n, 



-9(£ + k + lf 



(A4) 



Again n = \{£ — !){£ + 2), and in the start-up, S^ 2 = 
= e_2 and e_i =0. The origin value of the FDRK is 



(e+k) 



fe=0 



fc=0 



(A5) 

Despite the four-term recursion relation used to generate 
the series and its p-derivative, in using l|A5(l we have 
encountered no numerical instabilities. 



APPENDIX B: ASYMPTOTIC EXPANSION FOR 
OUTGOING SOLUTION 



To generate initial data for our numerical methods 
based on path integration in the complex z-plane (or 
sometimes the complex p-plane), we use the asymptotic 
expansion (|26J) about the irregular singular point at in- 
finity. Paper I describes this expansion for all the Regge- 
Wheeler cases, showing that the d^ w (a) obey a three- 
term recursion relation. For the Zerilli case, we write the 
expansion coefficients as df(a) — 5fe(cr)a' fe , then finding 
that the gk{c) obey the five-term recursion relation 



A k gk+i + B k g k+3 + C k g k+2 + D k g k+1 + E k g k = 0, (Bl) 



where 

A k = 

B k = 

C k — 

D k = 

E k = 



8crn 2 (fc + 4), 

An{k + 3)[6cr + n(k + 4)] - 8n 2 {n + 1) 

(k + 2)[18<7 - 4n 2 + (12n - An 2 ){k + 3)] - 12ti 2 

(k + 1)[(9 - I2n)(k + 2) - V2n] - 18n 

-9(k + l) 2 . (B2) 



fl k a k +3 + l k a k +2 + 5 k a k+1 + e k a k = 0, 



(A3) 



Here again n = \{t — 1)(£ + 2), and in the start-up 
C_3 = D_3 = = 0, D_2 = E-2 = 0, and £J_ a = 0. 
Typically, we have used fewer than ten terms in the ex- 
pansion (|26J) to generate initial data for numerical inte- 
gration, and have encountered no problems in using this 
expansion (despite potentially tricky issues associated 
with high order recursions). We have used dimensionless 
coordinates and Laplace rather than Fourier transform. 
Adjusting for these choices, I|B1|I and (|B2I) agree with an 
expansion given by Chandrasekhar and Detweiler in |34j . 
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